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In this paper, we study Bose-Hubbard models on the square and honeycomb lattices with complex 
hopping amplitudes, which are feasible by recent experiments of cold atomic gases in optical lattices. 

To clarify phase diagrams, we use an extended quantum Monte-Carlo simulations (eQMC). For the 
system on the square lattice, the complex hopping is realized by an artificial magnetic field. We 
found that vortex-solid states form for certain set of magnetic field, i.e., the magnetic field with 
the flux quanta per plaquette / = p/q, where p and q are co-prime natural numbers. For the 
system on the honeycomb lattice, we add the next-nearest neighbor complex hopping. The model 
is a bosonic analog of the Haldane-Hubbard model. By means of the eQMC, we study the model 
with both weak and strong on-site repulsions. Numerical study shows that the model has a rich 
phase diagram. We also found that in the system defined on the honeycomb lattice of the cylinder 
geometry, an interesting edge state appears. 

PACS numbers: 03.75.Hh, 67.85.Hj, 64.60.De 


I. INTRODUCTION 

Cold atoms in an optical lattice (OL) have been used 
as versatile quantum simulators for the last decade. In 
particular, the Bose gas system is highly controllable 
and has a rich phase diagram as a strongly-correlation 
system[lj. Recently, generation of an artificial mag¬ 
netic field in cold atom systems was experimentally suc¬ 
ceeded by using laser-assisted tunneling in a tilted optical 
potentialjl, Q, whose theoretical proposal was given by 
Jaksch and ZollerQ. These experimental methods can 
create a stronger magnetic field compared to that gener¬ 
ated by rotating optical lattice^, and therefore it is ex¬ 
pected that a strong-magnetic field regime corresponding 
to, e.g., the quantum Hall state is realized in cold atom 
systems. 

Bose gas system in a strong magnetic field has been 
studied very actively in the last several years. In par¬ 
ticular, study on two-component Bose gas system in a 
strong magnetic field in the continuum space predicts 
a bosonic analog of integer quantum Hall state (IQHS) 
for certain inter and intra-repulsions Generation of 
an incompressible vortex liquid was also suggested in a 
similar parameter region of Bose system in a strong mag¬ 
netic fieldQ. On the other hand for the Bose gas system 
on the optical lattice, the numerical exact diagnalization 
suggests an existence of a bosonic Laughlin state d B , 
and also a new kind of fractional quantum Hall state[lOj. 

There are many studies on the lattice boson systems 
with a dilute particle density and in a weak magnetic 
field[ni[i3. Formation of the vortex solid was predicted 
there. However for the system of interacting bosons in a 
strong magnetic field, the structure of the ground-state, 
in particular, structure of all the vortex-solid states are 
not known. Study on the complete global phase diagram 
of the system for various hopping amplitude, interactions, 
magnetic-field strength and particle filling, etc., is still 
missing. For lattice boson systems and certain related 


models in strong magnetic regime, a number of works 
have been reported so far. For classical spin models, Choi 
and Doniachfl^ found by an analytical method that the 
uniformly frustrated two-dimensional (2D) XY model has 
stable vortex-solid ground-states at / = 1/2, 1/3 and 
1/4, where / is the magnetic flux quanta per the funda¬ 
mental plaquette. Also, some of vortex-solid states were 
predicted by using Monte-Carlo simulations [l^ - [T^ and 
a Gross-Pitaevskii theory [TtL ; these studies showed 
some vortex-solid patterns for the 2D classical model at 
/ = 1/2, 1/3 and 1/4. Our previous study by using 
effective theory and Monte-Carlo simulations also pre¬ 
dicted some solid patterns at / = 1/2, 1/3 in two compo¬ 
nent lattice boson system assuming commensurate par¬ 
ticle density [l^. 

In theoretical models describing systems in a magnetic 
field, the hopping amplitudes acquire a nontrivial phase 
and become complex numbers. As a result, there appear 
various interesting phases, some of which are aforemen¬ 
tioned vortex solid, the bosonic Laughing state, etc. In 
the first half of the present paper, we clarify various vor¬ 
tex ground-states of a Bose Hubbard model (BHM) in a 
strong magnetic field by using extended quantum Monte- 
Carlo simulations (eQMC), in which effects of both phase 
and density fluctuations of the boson field are taken into 
account properly in the path-integral formalism. We 
found that a phase transition from the vortex solid to 
vortex liquid takes place as the on-site repulsion U is 
varied. In the second half of the paper, we study a 
bosonic analog of the Haldane model on the honeycomb 
lattice with the on-site repulsion U, which is sometimes 
called Haldane-Bose-Hubbard model (HBHM). Besides 
the nearest-neighbor (NN) hopping, the HBHM contains 
the next-NN (NNN) hopping with a nontrivial phase </>. 
The Haldane model was originally proposed as a fermion 
model that has the ground-state similar to the quantum 
Hall statefl^. In the recent paper[2lj, the HBHM with 
(/ = 7r/2 was studied by the dynamical mean-field theory 
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and an exact diagonalization of a small system. In this 
paper, we study the HBHM by the eQMC for both strong 
and weak on-site repulsion. We exhibit the global phase 
diagram, which should be compared with that obtained 
in Ref.[^, detailed critical behaviors of the phase tran¬ 
sitions, physical properties of an edge state in a cylinder 
geometry, etc. 

This paper is organized as follows. In Sec.II, we intro¬ 
duce the BHM in an artificial magnetic field and explain 
the derivation of the effective theory, which will be stud¬ 
ied by the eQMC. In Sec.Ill, we explain some details 
of the eQMC show the global phase diagram, which is 
obtained by the eQMC. In the phase diagram, there ex¬ 
ist various vortex-solid states and also vortex quantum 
liquid states for certain specific magnitude of the mag¬ 
netic field. Phase transition from the disordered state to 
vortex solid is studied by the finite size scaling and the 
critical exponents are estimated. In Sec.IV, the BHM in 
an magnetic field is studied by using a duality transfor¬ 
mation and origin of the vortex-solid patterns found in 
Sec.Ill is explained. Phase diagram obtained by vary¬ 
ing boson density is also shown. From these observation, 
possibility of bosonic Laughlin state is discussed by using 
a Chern-Simons theory. In Sec.V, the phase diagram of 
the HBHM is obtained, which has four phases. Detailed 
study of each phase is given by calculating the expecta¬ 
tion value of the current and phase correlation on links. 
Finally we investigate the HBHM in a cylinder geometry, 
in particular, we are interested in the edge state. Sec.VI 
is devoted for conclusion. 


II. BOSE HUBBARD MODEL IN A UNIFORM 
MAGNETIC FIELD AND THE EFFECTIVE 
MODEL 

In this section, we consider the BHM defined on a two 
dimensional (2D) square lattice. We start with the BH 
Hamiltonian in an artificial magnetic field with a vector 
potential A^{r), 

Hbyl =-J ^ -b h.c. 

+ Y.Unl ( 1 ) 

r 

where r denotes sites of the lattice and = x,y is the 
direction index and it sometimes denotes the unit vector. 
ar(aj.) is a bosonic annihilation (creation) operator and 
Ur = alar- J is hopping amplitude and U is on-site re¬ 
pulsive interaction. All these parameters are highly con¬ 
trollable in experiments [ij. A^(r) represents a uniform 
magnetic field perpendicular to the lattice plane. In the 
numerical calculation, we use 

where r = {x, y) and / is the magnitude of the magnetic 
flux per plaquette, and its range is 0 < / < 1 due to 


compactness of the phase degrees of freedom. The gauge 
field Afj_{r) creates the uniform magnetic field Bz = 
'^pAp = 27r/, where p denotes a directed close path 
around a plaquette and Ap is the vector potential on the 
path. This model is experimentally feasible^. 

In the following sections, we shall study the BHM in 
Eq.([T]) by the numerical MC simulations. To this end, we 
have to derive an effective model for the BHM with a pos¬ 
itive definite action. Although it is difficult to perform 
the direct quantum simulations because of the complex 
hopping in Eq. o, we can derive a useful effective model 
including relevant quantum effects by integrating out cer¬ 
tain degrees of freedom in the path-integral formalism. 

In previous work[T^, we derived the effective model 
for some related bosonic systems, and the effective model 
obtained there was numerically studied by the MC simu¬ 
lations. In this paper, we shall extend the previous meth¬ 
ods in order to search inhomogeneous states as a ground 
state. 

Let us start the derivation of the effective model men¬ 
tioned above. We first parameterize the boson variables 
in the path integral as = \/Pr + dpr where pr 

is the mean density at site r, which is regarded as vari¬ 
ational parameter in the MC simulation (see later dis¬ 
cussion). On the other hand, the variable 5pr represents 
quantum fluctuation of the density around the mean value 
Prj and then the Berry phase term in the action is given 
as {dr0{r))6pr, where the variable 0{r) is the phase of the 
boson field. By substituting the above parameterization, 
the partition function of the BHM defined by Eq. © is 
given as follows by using the imaginary-time r, 


ZBii= J[dSpr][de{r)]e-^^^, (3) 

•S'bh = J dr i{drO{r))Spr 

- y^^Jy'prPr-eu cos(6>(r) - 0{r + p) + Ap{r)) 
+jr{p, e)Sp + Y,{Upl + U5pl + 2Upr5pr)\ , (4) 

r ' 


where J'T{p,0)dp denotes the first-order contribution of 
the quantum fluctuation {Spr} in the hopping term J, 
and we have neglected the higher-order terms of {Spr} 
coming from the J-terms. As explain in the following 
section, we determine {pr} by the minimum-energy con¬ 
dition in the MC calculation. Therefore the terms of 
0{6p) in Eq.Q, except the Berry phase, cancel with each 
other. Then in Eq.Q, we perform the Gaussian integral 
for the density fluctuation Spr to derive the effective the¬ 
ory whose action is denotes by S^xy- Henceforth, we 
call this effective model Sqxy the quantum XY model 
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(qXYM), whose partition function is given by 

^qXY = J (5) 

5'qXY = f 

\ j. 

- J ^ ^PrPr+fi cos{9{r) - 0{r + p)+ Af,{r)) 

r,^ 

+ Y.Upl\. (6) 

r ' 

Here it should be remarked that S'qXY in Eq.® is real, 
and then the standard MC simulation is applicable for 
the numerical study of the system. Moreover it should 
be noticed that this model has the Lorentz symmetry and 
therefore it is expected that an excitation corresponding 
to the Higgs mode exists [l^, [l^. In fact this mode has 
been observed in optical lattice experiments . 

III. EXTENDED QUANTUM MONTE-CARLO 
SIMULATION 

In this section, we numerically study the qXYM by 
means of the eQMC explained in the previous section. 
In particular, we are interested in vortex dynamics in¬ 
duced by the strong artificial magnetic field with the 
vector potential A^{r), and obtain a global phase dia¬ 
gram of vortex states. For the study on phase diagram 
of the present system with the complex hopping ampli¬ 
tudes, the qXYM is quite useful as the phase degrees of 
freedom of the boson field plays an essentially important 
role. 

For the eQMC, we put the lattice spacing of the OL, 
ul, to the unit of length and also introduce a discretized 
lattice for the imaginary-time r with the lattice spacing 
At. Thus, the qXYM becomes a kind of 3D XY model 
defined on the space-time lattice, whereas its coefficients 
depend on the variational parameters {pr}- The lattice 
action of the qXYM is given as follows: 

•^qXY = E - ^)) 

r 

- JAry^^^prPr+f, cos(6>(r) - 0{r + p) + A^{r)) 

+ J2uatpI, 

r 

= 't)) +^^qXY, (7) 

r 

where r denotes sites on the 3D lattice. 

Here we briefly explain the eQMC to study the qXYM 
of Eq.([3). The qXYM action on path integral includes 
both the variational parameters {pr} and the dynami¬ 
cal phase variables 0(r)’s. We determines the variational 
variables {pr} by the minimum-energy condition by using 


MC methods. More precisely, we express the extended 
partition function [^qxv] hlq. ([7]) as 

[■^qXv] = J [dPr]ZqXYi{Pr}), 

^qXY({Pr}) = j [d0{r)]e-^^^'^. (8) 

In the practical calculation of Eq.®, we treat {pr} as 
slow variables in the MC local-update with keeping the 
mean value of {pr} constant. On the other hand, we 
perform the MC simulation for the dynamical variables 
0{r) as the ordinary variables for the fixed {pr}- From the 
experience, e.g. in Ref.jl^, we know that {pr} are quite 
stable under local updates for given values of parameters 
in the qXYM action. In the present study, we verified 
this stable behavior of {pr} for typical configurations of 
9{r), i.e., once {pr} is selected for fixed parameters of 
S'qXY by the MC method, {pr} does not change strongly 
for various configurations of 0{r). 

In the practical calculation, we employ the standard 
Metropolis algorithm with the local update [l^. The typ¬ 
ical sweep measurement is (50000—150000)x (5 samples), 
and the acceptance ratio is 40-50%. Errors are estimated 
from 20 samples by the jackknife method. Here we em¬ 
ploy the canonical ensemble, and therefore the mean to¬ 
tal number of boson, Pr^ is conserved during MC up¬ 
dates. This situation is suitable for real experiments of 
the OL[i|. We employ the periodic boundary condition 
and the lattice size is L^, where L is the linear system 
size. 

At has a relation, AtL = l/{kBT), thus its dimen¬ 
sion is 1/(energy). In our calculation, we set At = 1, 
regarded as unit of inverse of energy depending on sys¬ 
tem temperature. To identify the phase boundary, we 
calculate the internal energy E and “specific heat” C 
that are defined as follows, 

E = (i7qXY)/T", C = ((ifqXY - EL^f)lL\ (9) 

where (• • •) means the expectation value calculated as in 
Eq.®, 

(...)^ j[dpA{---){{Pr}). 

{■■■){{Pr}) = y'[d0(r)](---)e-^'-. (10) 

In FiglU we show the global phase diagram obtained 
by the eQMC for J = 5.0, and p = 2.0 (p is the mean 
density per site, i.e., p = {pr))- The phase diagram is 
shown in the (/ — U) plane and especially we focus on 
the magnetic flux regime 1/9 < / < 1/2. 

By calculating the Bose correlation function (a/or') by 
means of the eQMC, we have verified that for / = 0 and 
a sufficiently large hopping J, a SF with the long-range 
order of the phase 0{r) forms, i.e., {alar') —>■ c ^ 0 for 
|r —r'l —>■ 00 . As the value of / is increased, however, the 
SF order disappears at / ~ 0.02. This implies that the 
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/=l/9 



FIG. 1. (Color online) Phase diagram in f-U plane. The red- 
dotted lines indicates the f = Ijq vortex solid states and the 
blue-dotted lines the / = p/q vortex solid states. The phase 
boundary between the solid states and the vortex quantum- 
liquid states is determined by the behavior of vortex lines in 
the imaginary time direction. See FigO J = 5. 



FIG. 2. (Color online) The internal energy if as a function of 
the strength of the magnetic flux /. At / ~ 0.02, there exists 
a sharp peak that indicates a phase transition. Correlation 
function indicates that the system loses the SF at that point. 
J = 5 and U — 10. Calculated E indicates the existence 
of stable states (local minimums) for certain values of / like 
f = \ and whereas at f — ^ and a plateau appears. 

System size L = 12. 


critical magnetic field Be ^ 2 tt x 0.02. The internal en¬ 
ergy E shown in Fig l2] exhibits a sharp peak at / ~ 0.02, 
and therefore the phase transition seems to be of first or¬ 
der. As the value of / is increased furthermore, we found 
that at specific values of / like / = p/q, where p and q 
are co-prime integers, stable vortex-solid states form. We 
searched such a state for 1/9 < / < 1/2 and identified 
fourteen states. In the free electron lattice system in a 
uniform magnetic field, specific states for general f = p/q 


were predicted by Hofstadter[^. 

In Fig m the observed vortex solid states are indicated 
in the red- and blue-dotted lines for J = 5 and their 
snapshots are shown in FiglS) We define vorticity n(r') 
on the dual site r' of the site r as follows: 



sin(0(r -hi) — 9{r) + Ax{r)) 


-\- sin(0(r + X + y) — 9{r -I- i) -|- Ay{r -\- x)) 
E sin(6>(r E y) - 9{r E x + y) - A±{r -b y)) 


-b sin(6>(r) - 9{r + y) - Ay{r)) 


( 11 ) 


Expectation value of n(r') is calculated by the eQMC 
rather straightforwardly. 

We have found that the formation of vortex solid state 
depends on the lattice size L; i.e., for the f = p/q vor¬ 
tex solid state to form, the spacial lattice size must be 
qN X qM {N and M are natural numbers). As seen 
in Fig 131 the quantized vortices are pinned at sites of the 
dual lattice of the OL and they form solid pattern. In our 
previous study [l^ by the Monte-Carlo simulation assum¬ 
ing a homogeneous density of boson, some specific vortex 
solids (e.g., / = 1/2,1/3 and 1/4) were found. However 
in the present study using the eQMC, we found more 
general solid patterns / = p/q. This is due to the fact 
that the spatial density modulations are included as the 
variational parameter {pr} in the eQMC. As the den¬ 
sity fluctuations suppress the fluctuations of the phase 
degrees of freedom through the quantum uncertainty re¬ 
lation, and then the stable vortex solid formation is en¬ 
hanced in the present study. Behavior of the “specific 
heat” C as a function of the hopping J for / = ^ and | 
is shown in FigH) The results indicate the second-order 
phase transitions from the vortex solid to the disordered 
state as J is decreased. 

In Fig 131 we also show snapshots for incommensurate 
(‘irrational’) magnetic flux / ^ 0.23 and 0.27. It is ob¬ 
vious that vortices are located in sites of the dual lat¬ 
tice but they do not form a regular crystalline pattern. 
However they are rather stable against the MC updates. 
Therefore we conclude that an amorphous state of vortex 
forms at such fillings, although at the fillings belonging 
to the plateaus in Fig 131 a regular vortex lattice is to be 
maintained by the lattice pining effect. 

As we explained above, the density fluctuation influ¬ 
ences the dynamics of the phase degrees of freedom 9{r). 
Therefore it is interesting to see how the vortex behav¬ 
ior changes with the strength of the on-site repulsion U. 
Numerical simulation for the / = p/q vortex solid state 
shows that as the value of U increases, the vortex-lines 
in the imaginary-time direction starts to fluctuate. In 
FiglSl we show the typical vortex-lines in the / = 1/3 and 
1/4 cases. When U is small, U = 0.1, the vortex-lines 
are straight in the imaginary time direction, whereas for 
large U,U = 10, the vortex-lines are entangled with each 
other, i.e., the regular vortex lattices shown in Figl3]at 
some spatial layer tend to distort in adjacent layers and 



























5 





f^2/5(U=10) f=2/7(U=10) f^3/7{U=10) f^3/8(U=10) 


I =18 






f^2/9(U=10) f=4/9(U=10) f~0.23(U=10) f~0.27 (U=10) 


— L=12 
— L=18 
L=24 



FIG. 4. (Color online) Specific heat C for the system size 
L — 12,18 and 24. Scaling function of the FSS for / = ^ and 
f = - 

J 3- 


FIG. 3. (Color online) Snapshots of vortex fl(r') for p — 2.0 
and J = 5. They show vortex solid states in / = p/q (1/9 < 
/ < 1/2). The results indicates 14 patterns of the vortex 
solid states, which are all of possible combinations of coprime 
numbers p and q. For f ~ \ and 2, the rigid patterns appear, 
whereas for f = j and ^ locations of vortices are sightly loose. 
On the other hand for incommensurate /’s (0.23 and 0.27), 
vortices do not form a lattice. 

then the parallel-translated regular lattices reappear in 
some other layers. This different behavior of the vortex 
line stems from the quantum fluctuation, i.e., whether 
the degenerate vortex solid states in 2D are superposed 
or not through quantum tunneling processes. When the 
degenerate states in 2D are superposed, the state may be 
regard as a vortex quantum-liquid state. As shown in the 
phase diagram in FiglU the phase boundary between the 
vortex solid phase and the vortex quantum-liquid phase 
exists. However the phase boundary is not clear because 
the t7-dependent term in S'^xy in. Eq.Q generates only 
one-dimensional effect. Here we conclude that the on-site 
interaction U melts the rigid vortex-line states into the 
fluctuating vortex line states. In following section, the 
observed phenomenon in the above is studied by using a 
duality transformation of the qXYM. 

Finally, by using the finite-size scaling (FSS) of C, we 
numerically obtain the values of the critical exponents. 
See FiglU By the FSS hypothesis, C for the system size 
L, C'i(e), is parametrized as 

C'L(e) (12) 

where <l>(a:) is a scaling function and e = (J — Joo)/ Joo 
with the critical coupling for the system L —> oo. 
In Eg. (1121) . cr and v are critical exponent, in particular, 
v is the critical exponent of the correlation length ^ oc 


-^001 The results of the ESS are shown in FiglUfor 
f = 1I2 and 1/3. The critical exponents are estimated 
as 1 / = 0.80(0.55), cr = 0.44(0.24) for / = l/2(l/3) 

IV. DUALITY TRANSFORMATION 

In this section, we apply a duality transformation for 
the qXYM to understand the behavior of vortex observed 
in the previous section by the eQMC. This approach is an 
extension of the famous analysis on the 2D classical XY 
model, see for example Ref. [^ . 1^ . Duality-transformed 
model of the qXYM is described in term of the vortex 
density and topological current field. Besides the vortex- 
solid formation, it explains the vortex-line fluctuations in 
the T-direction controlled by the on-site repulsion U. Fur¬ 
thermore from the dual model, we find the / = 1/? rule 
for vortex solid pattern, which explains how the vortex- 
solid pattern is determined by the magnetic field /, the 
lattice size L, and vortex density-density interaction. 

A. Derivation of the dual-qXYM 

To derive the dual-qXYM, let us focus on the first and 
second terms in Eq. and define 

Zr+hop = J[d9{r)]e-^-+^'>-, (13) 

Sr+hop = - XI cos(0(r + t)- 9{r)) 

r,r 

- X cos(9(r -hp) - 9(r) -h A^(r)), (14) 

r,/j, 
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V 


X 


/ = l/3, t/ = 0.1 


/ = l/3, C/ = 10 






/ = 1/4, U = 0.1 


/ = 1/4, f/ = 10 


FIG. 5. (Color online) The typical vortex lines in / = 1/3 and 
1/4 along the imaginary time direction r. The red dots are the 
vortices residing in the dual site r. As the on-site repulsive 
interaction U increases, the vortex-lines are entangled, and 
the quantum tunneling from the vortex lattice shown in FigO 
to another parallel-translated one takes place. When U is 
small, the vortex solid states does not change along imaginary 
time direction. 


where we have set At = 1, = (pr) = 1 for simplicity. 

(Effects of fluctuations in the local density Pr will be 
discussed in Sec.IV.C and Sec.IV.D.) In the following, 
we express the action Sr+hop as follows for notational 
simplicity, 

Sr+hop = - cos(0(r + p) - 9{r) + A^{r)), (15) 

where we have redefined direction labels, the couplings, 
and vector potentials as 

/i = 0, l(x),2(y), Jo = ^, Ji = J 2 = J (16) 
and use the gauge 

Ao{r) = 0, Ai{r) = 2Trfy, A 2 {r) = 0. (17) 


Let us apply the periodic Uaussian approximation 


Eq.jni) [13,11 


cos(e(r+p)-e(r)+Aij,{r)) 


00 

E 

nu. — — oo 


exp 


Ja 


{9{r + p) - 9{r) 


+A^{r) - 27rn^(r')^ 


(18) 


where the vector field ^^(t') is defined on the sites r' of 
the 3D dual lattice, and this field represents 27r period¬ 
icity of the cosine term in Sr+hop- Substituting Eq. (fT51) 
into the partition function Eq. (1131) , we obtain the follow¬ 
ing equation by using Poisson formula (see for example, 

Ref. dal), 


-'T-\-hop OC 


OC 


f [d9{r)] X] exp - ^ 

1 L r.u 


G(r') ■- ’■-M 

;2 (^1 


2Ju 


- ilp{r'){9{r + p) - 9{r) -h A^{r)) 


E E “ ilpir)Apir)] n^i5(^ - l^{r' - p)), 

G(r') ^ 


(19) 


where we have used the Hubbard-Stratonovich transfer- gauge fields d^(r), with which lp{r') is expressed as 

mation, and the vector integer field lp{r') has been in- ^+{ 1 '') = vd-x{r'). By substituting the above ex- 

troduced as the boson-current variables. To solve the 5- pression of lp{r') into the action, we have 

function constraint in Eq. (1191) , we introduce new integer- 


Zr+hop OC ^ ) OXp f ^ ) 

5,^u(r) ^ r,p 


{ipuxy udx{r)Y 
AttJu 


- KpuxVvdx{r)A^{r) 


( 20 ) 


where the operator is lattice nabla, defined by 


^pA,j{r) = Ay{r + p) — Ay{r). Here we use again Poisson 
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summation formula, and then the integer variables a^{r') 


J 


are transformed into the continuum variables and 

the partition function is expressed as, 


Zr+hopO^ / [#^(r')] ^ expl-^ 




AttJu. 


- v(l)x{r')A^{r) + 27rzTO^(r')(/)^(r') 


/ Wpir')] X! - 2mf^4>fj,{r') + 27rzw^(r')(^^(r') 


m^(r 

/ oo 

-OO 


( 21 ) 


with 


Sv{jp{r'),(i3^{r')) = ^ 
r' ,/j. ^ 


1 

47rJ^ 




( 22 ) 


where we have introduced the three-component flux 
field fp{r) for the notational simplicity 27r/o(r) = 

eoijWiAj{r), /i( 2 )(r’) = 0, and j^(r') = 27r(TO^(r') - 
fp{r)). By integrating the “gauge fields” we ob¬ 

tain the final expression of the dual model. 

As the action 5'„ in Ea. (E^ is invariant under a local 
gauge transformation (jxxir') —^ 4>xir') — V\a{r') with 
an arbitrary scalar function a{r'), we have to fix the 
gauge of (/i'(r') to integrate out To this end, we 

consider the case of = J for simplicity and employ the 
Lorentz gauge 1 1 = 0. In the present path- 

integral formalism, this gauge condition is easily imposed 
by adding the gauge-fixing term (^^ V^(()^)^ to the ac¬ 
tion S^. In the Lorentz gauge, 

Sv - + (j - terms) 

r' 

= + (■? “ terms). (23) 

r' ,ii,X 

Integration over (f>{r) in Eg. (1211) can be performed 
straightforwardly and obtain the dual model of the vor¬ 
tex density and vortex current j^(r') as follows. 



S'duai = -jYl \r,r')jf,{r'). (24) 

r',r" fi 

In what follows, we shall study the stationary state of 
the vortex configurations like the vortex solid. We first 
interested in the vortex density-density interaction in the 
dual model, i.e., the interaction between the component 
jo in Ea. (|24|) . and evaluate the corresponding term as 

- 7rJ^jo(r) ^ jo(rO = nJ jo{k)-Yjo{-k), 

( 25 ) 


where jp{r) = z (B.Z. refers to the first 

Brillouin Zone) and = -^ sin , xhe same re¬ 
sult with Ea. (l25|) is obtained in the Coulomb gauge 
X)a=i 2 = 0. In the long-distance region |r'—r"| 

ol, Eg. (051) behaves as 

'^J^jo{k)Yjo[-k) 

B.Z. 1^1 

r',r" k 

oc ttJ ^ joir')log\r' - r''\joir") +a^^io(?'')^ 

= ttJ ^ ijoir') - 27r/) log|r' - r"|(jo(r') - 27r/) 

+ «(^^Oo(»'') - 27r/)^ . (26) 

where jp,{r') = 27rm^(r'), and we have introduced the 
a-term by regularizing the infrared divergence in the k- 
integration and then a ^ ln(L/aL). It is easily veri¬ 
fied that the current-current correlations (ji — ji) and 
(j 2 — j 2 )-terms have a similar form with Eo. (l26E but its 
coefficient is l/4t7. Therefore as U is getting large, the 
current correlation becomes weak and nontrivial config¬ 
urations of the current j = (ji,^ 2 ) appear. This means 
the correlation of vortex density in the r-direction is get¬ 
ting weak as we observed in the numerical simulations 
in FiglHl In following section, we explain the generated 
patterns of the vortex solid by the eQMC from Ea. (l26l) . 

B. Vortex solid / = 1/? rule 

In previous section, we derived the dual-qXYM model. 
For the stationary configurations, the vortex-density in- 











teraction has the logarithmic form similar to that of 2D 
XY model, and its coefficient is proportional to the hop¬ 
ping amplitude J. Moreover, the infrared singularity 
gives the “charge-neutrality condition” as 

oc (^(27rTOo(r) - 2TTf))^ 0. (27) 

r' r 

Thus for sufficiently large hopping J, the BEC is real¬ 
ized in the system and the above condition decides the 
total number of vortex in the system (more precisely, the 
number of (vortex — anti-vortex)). Here, we again show 
the vortex-density interactions, 

J ^(27rTOo(r) — 27r/) log \r — r'|(27rmo(7’') — 27r/) 

r,r' 

+ a(^{2Trmo{r)-2Trf)f. (28) 

r 

From the above interactions, we qualitatively explain the 
patterns of vortex-solid observed by eQMC in the previ¬ 
ous section. For simplicity, we shall consider the case of 
the magnetic flux f = 1/q. 

1. In the neutrality condition Ea. (l?7l) . / = 1/q and 
mo(r)/27r is an integer, q adjacent plaquettes have 
to be grouped for making qxl/q = 1 flux quantum. 
To achieve this globally, the lattice system size has 
to be qNx X qNy {Nx and Ny are natural numbers). 
This is the lattice size condition of the system. 

2. The lattice system is divided into a number of q- 
adjacent plaquettes. One quantized vortex resides 
in arbitrary one of g-adjacent plaquettes. To re¬ 
alize low-energy configurations, the vortex density 
interactions have to be minimized. This decides the 
distribution pattern of a number of vortices. 

3. The number of the distribution pattern corresponds 
to the number of the degenerate state of the / = 
1/q vortex-solid state. When the particle density 
is high enough, a superposed vortex-solid patterns 
is realized because of the long-range repulsion be¬ 
tween vortices. 

The above consideration can be easily extended to the 
more general cases with f = p/q {p and q are co-prime), 
and from this consideration it is expected that all possi¬ 
ble vortex solids can be observed by the present eQMC 
simulations 
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FIG. 6. Numerical results of the dilute density regime. (Left 
panel) The phase diagram in the p — f plane. The red bro¬ 
ken lines show the observed vortex solid states in which the 
vortices are pinned on the dual lattice. (Right panel) Snap¬ 
shots for / = 1/3. At p = 2, the robust vortex solid forms. 
As the density is decreased, ‘melting’ of solid takes place. At 
present, we do not have clear understanding of the states in 
the dilute regime. However, we expect that they are a kind of 
vortex-liquid state. In the right column, we show the results 
for / = 1/2 and p = 0.25, in which a bosonic fractional Hall 
state is expected to form. The state is rather homogeneous 
but slightly unstable against the MC updates. 


disappear. As shown in FigjSJ locations of vortices do not 
exhibit any spatial patterns in the low-density region and 
the transition may be interpreted as a melt of the vortex 
solid to a vortex liquid Q- Specific heat C shows a sharp 
peak at the transition point for some values of /, but for 
other values of / there is only a moderate peak in C. 
Due to the diluteness of bosons, the effect of the lattice 
is weekend and the system approaches to the continuum 
system. Pining of vortex by the lattice is less effective, 
and the topological vortex number is spread rather wide 
spatial regions. These findings are in agreement with the 
intuitive picture such that for the low-density limit, the 
lattice works simply as a mesh employed for numerical 
study. In fact, the low-density limit means small hop¬ 
ping parameter J{prPr+y)^^‘^ and its relatively large lo¬ 
cal fluctuations, and then the interaction term in Ea. ((28ll 
becomes less effective. This observation is quite instruc¬ 
tive for the discussion on the realization of the bosonic 
Laughlin state in the following section. 


C. Phase diagram in dilute boson regime 

In Figini we show the phase diagram of vortex state in 
the low density regime. As explained in Ref. [2^, recent 
experiments can decrease the atomic density without los¬ 
ing the phase coherence by applying a micro-wave to local 
regimes of the OL and blowing away excited atoms from 
the OL. By eQMC calculation, we found that as the aver¬ 
age density of bosons is decreased, the vortex-solid states 


D. Bosonic Langhlin state and Chern-Simons 
theory of the Bose Hubbard Model 

In recent years, it was suggested that a state similar to 
the bosonic Laughlin state forms in the Bose gas system 
loaded on the OL. In the first-quantization language and 
also in the continuum space, the bosonic Laughlin state 
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is given by the following wave function [s^, 

■ ■ ,Zn) = (29) 

j<k 

where Zj = Xj + iyj and m is an even positive inte¬ 
ger. From Eg. d^^ . it is obvious that bosons in the 
Laughlin state have the hard-core nature as —?> 0 for 
Zj —>■ Zk (j ^ k). In the previous work[^ Q, a numer¬ 
ical diagonalization of a small system showed that the 
ground-state of the hard-core BHM has a good overlap 
with the wave function in Eg. (1291) for a weak magnetic 
field / < 0.2 and low density p < 0.1 case. 

It is interesting to ask if a state that has similar proper¬ 
ties to the Laughlin state forms in the boson gas system 
on the OL for relatively large /’s and particle density. 
Such a state has a certain hidden order and low-energy 
excitations behave like anyons. To study this problem, 
Chern-Simons (CS) theory for the fractional quantum 
Hall effect (EQHE) is quite useful. We previously in¬ 
troduced a lattice version of the CS theory that is well 
suited for studying the present system[3l|. 

By using this formalism, we transform the original bo¬ 
son operator Ur in Eq. © to another boson operator br, 
which we call CS boson, by attaching m-magnetic flux to 

dr , 


dr — ? 

Vr = exp 


im 


'^6{r,x'){al,ar>) , 


(30) 


where x {x') denotes a site of the dual lattice paired 
to site r {r') of the original lattice, and 9{r,x') is the 
azimuthal angle function on the lattice. Please notice 
ojor = blbr- Then, let us consider the specific case 
/ = mp where p is again the average particle density. In 
this case in terms of the CS boson br, the BHM Hamil¬ 
tonian in Eq. is expressed as, 


Hbk =-J ^ blw}Wr+f,br+^ + h.c. 

+ '£U{blbr)\ 

r 

ITr =exp d(r, a;')(6|!,6r'— p) , 

x' 


(31) 

(32) 


where we have employed the symmetric gauge for A^{r), 
and used the identity 27re^yVyG'(x, x') = V^d(r, x') 
[ei 2 = —£21 = 1, £11 = £22 = 0] with the two-dimensional 
lattice Green function G(x,x'). Beautiful CS gauge the¬ 
ory can be constructed for the system Eqs. m and dsa) 
in the Lagrangian formalism, but here we only discuss the 
possible mean-field (ME) solution to the ground-state of 
the above system. 

The ME analysis indicates the naive candidate for the 
ground-state of Hbh in Eo. d^ such as b^^ = p^^^, i.e., 
in which the BEC of the CS boson br takes place. Erom 
Eg. (1501) and f = mp, this MF solution accompanies the 
homogeneous ‘condensation’ of the 27r/ flux quanta per 


plaquette to cancel out the external magnetic field. This 
MF state corresponds to the bosonic Laughlin state and 
therefore in the bosonic Laughlin state, the grand-state 
has the homogeneous particle density p and also the ho¬ 
mogeneous vortex density 27r/ per site and per plaquette, 
respectively. However for sufficiently large p and /, the 
above condition is not satisfied by the configuration of 
ar because the duality transformation shows that vortex 
density at each plaquette TOo(r') takes an integer value. 
This is nothing but the lattice pining effect of the vortex. 

On the other hand for small p and /, the fluctua¬ 
tions of the parameter J{prPr+fj,y^‘^ cannot be neglected 
and therefore the direct application of the result of the 
duality-transformation is questioned. It is quite useful to 
study, for example, the case with / = 0.5, p = 0.25 and 
therefore m = 2. Snapshot of the topological number 
(vortex density) is shown in FiglO] It is obvious that a 
genuine homogeneous configuration is not realized there 
but topological number tends to smear compared with 
the higher-density case in Figsl3]and[ni This result seems 
to indicate the possibility of the bosonic Laughlin state 
for very low particle density as indicated by the work 
on the very small systems in Ref[^, 0 . Methods in the 
present paper are also applicable for the study of such 
low density cases, i.e., the hard-core BH model. To this 
end, the slave-particle representation for the hard-core 
boson is quite useful This problem is under study 
and we hope that the results will be reported in another 
publication in future. 


V. HALDANE BOSE-HUBBARD MODEL 
STUDIED BY eQMC 

Resent experiments on clod atoms succeeded in real¬ 
izing the Haldane-Bose-Hubbard model in a honeycomb 
optical latticed^. This system is a bosonic analog of 
the celebrated Haldane mo del [13 , which is a fermionic 
system and possesses a topological phase due complex 
hopping amplitudes on the honeycomb lattice. Besides 
the NN and NNN hopping terms, the on-site repulsion 
exists in the HBHM, and as result of the competition of 
these three terms, the model has a rich phase diagram. 
The HBHM realized by cold atomic gases on the opti¬ 
cal lattice is highly controllable, e.g., besides the average 
particle density, the on-site repulsion, the hopping am¬ 
plitude and the artificial magnetic flux generated by the 
NN hopping can be controlled. 

In this section, we shall study the HBHM by means of 
the eQMC and obtain the global phase diagram. To this 
end, we calculate the internal energy, specific heat and 
certain correlation functions to identify order of phase 
transitions and physical properties of phases. By using 
the FSS, we obtain the critical exponents for phase transi¬ 
tions. The previous work[2l| clarified the phase diagram 
of the HBHM, however the system size in the numeri¬ 
cal calculation (the exact diagonalization) is small, and 
only the case of unit filling was considered. Therefore 
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the present study using the eQMC is complementary to 
the previous work. In particular, besides those found in 
Ref.(2l| we have found another phase boundary. 

As in Ref. [2l| , the Hamiltonian of the HBHM is given 
by 

Hhbh = -^1 ~ J-2 e“*‘^a|aj 

{i,3) (dd» 

+U'^a\aia\ai, (33) 

i 

where (j,j) denotes NN sites, {{i,j)) NNN sites of the 
honeycomb lattice. The parameters Ji and J 2 are the 
NN and NNN hopping amplitudes, respectively, and (f) 
is a constant phase of the NNN hopping. By a calcula¬ 
tion similar to that in Sec.II, the effective action for the 
HBHM is obtained as 

5'hbh = J 

^ i 

- Ji y] ^ypiPj cos{e^ - 9j) 

(i,i) 

- J 2 ^ Vft^cos(6'j -6j+(j))\. (34) 

{(.ip)) ^ 



FIG. 7. (Color online) Honeycomb lattice with the periodic 
boundary condition (torus) and also that of the cylinder ge¬ 
ometry. The system size can be defined from A sub-lattice 
size, its lattice is a square L x L lattice. The imaginary-time 
T lattice is fixed Lr = 8 and has periodic boundary condition. 

For the numerical calculation, we consider the follow¬ 
ing system; 

1. We consider the unit-filling case and in most of 
the calculations we set the variational parameters 
Pi = p = 1 for simplicity. Full MC simulation of 
the action (IM)) verifies the validity of the above as¬ 
sumption. 

2. We put (j) = '^/2 as in Ref.[2l|. 

3. Both the torus and cylinder geometries are consid¬ 
ered. See FiglT] 



FIG. 8. Phase diagrams of HBHM with U = 0.1 (top left) 
and 1/ = 10 (top middle). Galculations of the various physical 
quantities are shown, which are used for identihcation of the 
phases. Between the SF and CSF, there exists a coexisting 
SF-I-CSF phase. “Specific heat ” C along the line (c) in the 
phase diagram (bottom right) exhibits two peaks at Ji ~ 1.6 
and 2.3. The system size is L = 6 


4. To classify the physical meaning of the observed 
states, we measure the current operators defined as 

A, = 2Im(j,,(e*(®*-''^))), (35) 

and also the link correlation defined by, 

i?,, =2Re(j,,(e*(®--®^))), (36) 

where Jij is Ji and J 2 e*‘^ for the links connecting 
NN and NNN sites, respectively. 

As explained in the previous work Ref[2l|, for small Ji/U 
and J 2 /U, the system is the plaquette Mott insulator 
(PMI) that has only finite local correlation ^ 0 

and vanishing nonlocal correlations. On the other hand 
for a large Ji J 2 , the ordinary superfluid (SF) forms 
and it has a positive expectation value > 0 and 

= 0. For the weakly interacting limit, “ 

2pJ2. For Ji <C J 2 , the J 2 -term in the Hamiltonian 
Eg. (1551) dominates and the chiral SF (CSF) forms with a 
negative expectation value < 0. For the weakly 

interacting limit, ~ —‘^pJ 2 - 

In this work, the HBHM is studied by the eQMC and 
the phase boundaries are identified by calculating the in¬ 
ternal energy E and the “specific heat” C as in the study 
on the BH model in the previous section. We furthermore 
calculate the density of states, N (S'), to identify the order 
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J X 

'’2 

(a) PMI-CSF 



FIG. 9. “Specific heat” C and the scaling function of the FSS 
<l?(a;). Upper panels, the PMI-CSF phase transition. Lower 
panels, the PMl-SF phase transition. FSS indicates that both 
phase transitions are of second order. U = 10. 


of phase transitions, which is defined as 

^HBH = J [dd(r)]e-®«^« 


= dS e 


-s 


[d0{r)] S{S — Smbh) 


= J dS e-^N{S). 


(37) 


are induced by growing density fluctuations. “Specific 
heat” C and the scaling function of the FSS are shown 
in Figl3 The results indicate that both the PMI-SF and 
PMI-CSF transitions are of second order. The critical ex¬ 
ponents are estimated as v = 0.95(1.25), cr = 0.42(0.26) 
for PMI-SF (PMI-CSF). On the other hand for the phase 
transition from the SF to the SF-I-CSF phase and also 
from the CSF to the SF-I-CSF phase, the calculated C 
for various system sizes does not exhibit the FSS, nor 
the density of states N{S) defined by Eg. (1571) exhibits a 
double-peak shape at the phase transition point. More 
detailed study using a large-scale systems is needed to 
clarify the properties of the phase transition SF (CSF)—^ 
SF-bCSF. 



SF+CSF CSF 

Ji = 2.2 Ji = 0.6 

J 2 = 1.4 J 2 = 1.4 


1.0 




-0.5 


FIG. 10. Snapshots of vortex Ha and Hs in the SF, CSF 
and SF-I-CSF phases. In the SF, phase of BEC has a uniform 
distribution. On the other hand in the CSF, 120°-structure 
forms. In the SF-I-CSF phase, a spatial mixture of them ap¬ 
pears. 


On a phase transition point, we calculate N{S) by the 
MC simulations. If e~^N(S) has a double-peak shape 
as a function of S, the phase transition is of first or¬ 
der, whereas a single-peak shape of e~^N{S) indicates a 
second-order transition. 

We first consider the system with the periodic bound¬ 
ary condition, i.e., a torus. FigjS] exhibits the global 
phase diagrams of the HBHM with U = 0.1 and 10. In 
the phase diagram, there are four phases, three of which 
are found in the previous work 2l|. i.e., the PMI, SF, 
and CSF. The fourth phase is a coexisting phase pos¬ 
sessing both the SF and CSF correlations. Calculated 
“specific heat” C is shown in Figl 8 ]as a function of Ji, 
which indicates the existence of two phase transitions for 
J 2 = 1.2 and U = 10. The NN and NNN currents are cal¬ 
culated for the identification of the phase and the results 
are shown in FiglU In the SF, the intra-sublattice cur¬ 
rent > 0, whereas in the CSF < 0 as explained 

in Ref. . In the SF-I-CSF coexisting phase, the value 
of I A A, Iaa^ is IaT < dAA_ < Phase diagram of 

other values of U has qualitatively the same structure 
with that in Fig 15] However, as decreasing on-site inter¬ 
action U, the phase boundary shifts because of the sup¬ 
pression of quantum fluctuations of the phase di, which 


It is interesting to see if a vortex solid forms in the SF 
and CSF phases, and in particular how the vortex solid 
is deformed by the coexistence of the SF and CSF in the 
SF-I-CSF phase. To this end, we calculate the triangle 
vorticity of the A and H-sublattices, HA(B)(r), which is 
defined by 

(?■) = ^ ( sin( 6»2 -01 +^) + sin( 6'3 - 02 + (/)) 

+ sin( 6 (i - 03 - 1 - ^)'j , (38) 

where the detailed definition of 0’s, see Fig [TUI and (j) = 
7 r /2 in the present case. For the configuration of the 
uniform 0i{r), Ha(b)(^) = whereas for the 120 °- 
configuration, flA{B)(j’) = ~l/2. Snapshots of flA(B)(r') 
are shown in FigfTUl and the results indicate that the 
uniform-0 configurations are realized in the SF, whereas 
nearly 120°-configurations are in the CSF, as expected. 
On the other hand, the SF-I-CSF phase has a nontrivial 
distribution of the vortex. From this result, we expect 
that immiscible SF droplets exist in the CSF and vice 
versa in the SF-I-CSF state.. 

Finally, let us study the HBHM in the cylinder geom¬ 
etry shown in Fig[71 In particular we are interested in 
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boundary 
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FIG. 11. Phase diagrams of HBHM in the cylinder geome¬ 
try. Calculations of the various physical quantities are shown, 
which are used for identification of the phases. Near zigzag 
edges, the quasi-one-dimensional SF state forms in the bulk 
CSF state. Upper panel for U = 0.5 and lower panel U — 10. 
The system size is 1/ = 6 


BEG behavior near the zigzag edges of the cylinder. In 
recent experiments, such a sharp edge boundary can be 
created by using the optical-box trap method [3^. For 
the SF phase of the HBHM, Bogoliubov excitations near 
the edges was recently studied[35l|. On the other hand in 
the present study, we focus on the SF-hCSF phase that 
was observed by the eQMC in the previous section. It 
is expected that near the edges a quasi-one-dimensional 
excitation forms and the Ji-term in Eg. (1331) dominates 
there. Therefore, the SF state appears near the edges of 
the cylinder. This expectation is verified by calculating 
lij and Rij. See FigfTTl Schematic picture is that for the 
bulk CSF state, the SF state appears near the edge of the 


cylinder BEG, whose bnlk state is the CSF for J 2 > J 2 C, 
where J 2 C is the critical value of J 2 . On the other hand 
for smaller J 2 < J 2 C, the whole system is the SF. 


VI. CONCLUSION 


In this paper, we studied the BHM and HBHM with 
complex hopping amplitudes that are recently realized 
by experiments of the cold atoms [3, H, We first ex¬ 
plained the numerical methods that we call the eQMC. 
Most of the numerical results in the present paper have 
been obtained by the eQMC. 

For the BHM in various magnetic fields, we focused on 
the vortex-solid state and clarified its existence for vari¬ 
ous flux quanta / = p/q. Various spatial patterns of the 
vortex-solid have been identified by calculating the wind¬ 
ing number for each lattice plaquette. We furthermore 
found that the on-site repulsion U plays an important 
role for the phase transition from the vortex solid to vor¬ 
tex liquid. By the duality transformation of the effective 
mode of the BHM, we discussed the condition on / for 
the formation of the vortex solid. 

For the HBHM, we have obtained the global phase dia¬ 
gram of the ground-state with the complex NNN hopping 
with (j) = 7r/2. Besides those found in Ref.[2l[, we have 
found another phase in which the SF and CSF regions co¬ 
exist. Uniform and 120° structures of the phase degrees 
of the freedom of the condensed boson field are ‘entan¬ 
gled’ with each other. More detailed study is needed to 
clarify the dynamical properties of the SF-I-CSF phase 
like a quantum (in)stability, etc. This is a future prob¬ 
lem. Finally, we studied the SF-I-CSF state of the HBHM 
in the cylinder geometry, and found that the quasi-lD SF 
state appears near the edges of the cylinder. This phe¬ 
nomenon might have some connection to the edge state 
in the quantum Hall like state. This is also a future prob¬ 
lem. 
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